A statistical approximation to solve ordinary differential equations 
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We propose a physical analogy between finding the solution of an ordinary differential equation 
(ODE) and a N particle problem in statistical mechanics. It uses the fact that the solution of 
an ODE is equivalent to obtain the minimum of a functional. Then, we link these two notions, 
proposing this functional to be the interaction potential energy or thermodynamic potential of an 
equivalent particle problem. Therefore, solving this statistical mechanics problem amounts to solve 
the ODE. If only one solution exists, our method provides the unique solution of the ODE. In case 
we treat an eigenvalue equation, where infinite solutions exist, we obtain the absolute minimum of 
the corresponding functional or fundamental mode. As a result, it is possible to establish a general 
relationship between statistical mechanics and ODEs which allows not only to solve them from a 
physical perspective but also to obtain all relevant thermodynamical equilibrium variables of that 
particle system related to the differential equation. 
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I. INTRODUCTION 

Ordinary differential equations are one way physicists 
model physical phenomena. To solve them, a great vari- 
ety of analytical and numerical methods are available in 
the literatureii^. 

Trough this article we wish to solve some problems 
that arise in vibration theory making an analogy of them 
with a system of particles, a well known issue treated in 
Statistical Mechanics. If this analogy is made, all the 
tools used in Statistical Mechanics are available for tack- 
ling these type of problems. For example, one is able to 
calculate the entropy or the specific heat of the system, 
something which is, at least, difficult to imagine in a typ- 
ical vibration problem. From a more general perspective, 
we extend the same analogy and apply it to find the solu- 
tions of ordinary diff'erential equations (ODEs) of general 
type (see section HI)) . 

The general method proposed here could be resumed 
in the following way. From a mathematical viewpoint, 
it is possible to obtain the solution of an ODE from 
the minimization of a functional. On the other hand, a 
thermodynamic potential is represented by a functional. 
We connect these two notions, proposing the functional 
which leads to the solution of an ODE, to be the interac- 
tion potential energy or thermodynamic potential of an 
equivalent particle problem. Therefore, solving this sta- 
tistical mechanics problem amounts to solve the ODE. 
For the case where only one solution exists, the method 
allows to obtain the unique solution of the ODE. Clearly, 
when we treat vibrational problems, an eigenvalue equa- 
tion must be solved with infinite possible solutions. This 
case our method is capable of obtaining the fundamental 
mode or absolute minimum of the corresponding func- 
tional. 

The principal aim of our work is to show (a) how the 
previously proposed correspondence is used to attain the 
solution of ordinary differential equations, (b) the theo- 
retical basis upon which it is sustained and (c) the way 
it can be treated within a stochastic framework. 



The present work is organized as follows. In section 
im we show how the problem of solving a general kind of 
ODE can be transformed to a problem of minimization of 
a functional (weak solution of the differential equation''). 
Section IIIII presents the foundations of equilibrium sta- 
tistical analysis viewed as a minimization of a thermo- 
dynamic potential. This provides the theoretical the ba- 
sis of our developed method. Section IIVI describes the 
general aspects of the numerical scheme and algorithm 
implementation. Section |V] is devoted to find the solu- 
tion of an ODE employing a stochastic framework. This 
method allow us to solve the ODE proposing an equiva- 
lent Langevin (stochastic) equation. At the end, in sec- 
tion lVIl several numerical experiments which will sustain 
our theoretical hypothesis will be developed and solved. 
Finally, some conclusions will be outlined. 



II. MATHEMATICAL PROBLEM 

There are many possible approaches to find the solu- 
tion of ODEsii^. From all of them, we are interested in 
that one which transforms the problem of finding the so- 
lution of an ODE into the minimization of a functional. 
To do this we must recall some results of functional anal- 
ysis. 



A. Obtaining a minimum principle for the 
differential equation 

The most general form of an ODE of order 2k can be 
written as^ 



^«(:r) = ^(-l)Xp.?/(^'(^))^^' -/ (1) 

1=0 

where u{x) is an element of a Hilbert space "Ha, a; G 5R 
and u^'^\x) is the derivative of order (z) of u{x). Also 
/ e L2{a,b) {L2{G) Hilbert space of functions square 
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integrable in the domain G) and Pi{x), i = 0, 1, are 
functions continuous with their derivatives up to the i—th 
order inclusive in the closed interval [a, b] which satisfy 

p^{x) > in[a, b], i = 0, 1, k-1 (2) 

Pk{x) > p in[a, b], p = const. (3) 

For the rest of the paper, homogeneous boundary condi- 
tions will be considered with no loss of generality. 

Then, it is possible to look for the generalized solution 
u* of equation ([1]) as that element which minimizes the 
functional 




in the corresponding Hilbert space Ti.A, where oper- 
ator A is defined. There exists a requirement of posi- 
tive definiteness of operator A that ensures that the ex- 
treme value of (|3]) or (O is the desired solution u* of the 
probleroii. 

For the eigenvalue problem, Au ~ Xu — 0, A e SR, 
the functional to be minimized coincides with the corre- 
sponding eigenvalue A and it reduces to 

F{u) ^ A = ^ ° (5) 

/ u-^dx 
J 

In this case, infinite minima exist, each of them corre- 
sponding to the associated eigenvalues 

Ai < A2 < ... (6) 

Equation ^ states that the first (or fundamental) eigen- 
value Ai is the minimum possible eigenvalue and the cor- 
responding eigenfunction is ui; i. e. F{ui) is the absolute 
minimum of F{u)^. 

Summing up, a differential operator can be trans- 
formed into a functional so that minimizing this func- 
tional equals to solving the ODE. 

III. STATISTICAL TREATMENT 

As a second step towards attaining the solution of 
ODEs, we must adequate the problem to be treated as 
a standard particle problem in statistical mechanics. To 
this end, we must think on a n particle system in a one 
dimensional box of size L, where a point in the n di- 
mensional configurational phase space can be denoted 
by q = {^ii 92, ■•■<?«}, being qj the position of particle 
j. Now, a link between the representation of this system 
in configurational phase space q with the description of 
possible candidate solutions of the ODE u{x) (see equa- 
tion [1]) is proposed through an adequate discretization of 
u{x). 



A. Discretization of u{x) 

Consider now any possible candidate to the solution of 
a particular kind of ODE, named u{x). This continuous 
function is an element in the corresponding Hilbert 
space Ti-A-, a; G Sft where operator A is defined, and 
belongs to the domain of the functional F{u). The 
proposed discretization of u will be settled from the 
consideration of n sites labelled with the discrete index 
Z of a 1 dimensional lattice. Then, the function u(x) is 
taken as the limit n — > 00 from the field functions ui{xi). 
So the link between our system of particles and that 
element u{x) is naturally settled as ui{xi) = qi. Each 
ui can take any value between —00 and 00. A possible 
configuration of our n particle system is then considered 
by specifying the value of ui at each site I or simply 
the lattice vector q = {wi(a;i), u/(x/), ...u„(x„)}. 
Through this connection, it is possible to apply the 
whole statistical framework to treat this equivalent new 
problem. 



B. Equilibrium statistical analysis 

The next step is to formulate the problem as in stan- 
dard statistical mechanics. To do this, we will consider 
our particle system in the canonical ensemble. So far 
we have made a connection between the domain of the 
functional F{u) with the microstates of a thermal sys- 
tem, represented by q. With this in mind, and using the 
expression for the probability density in the canonical 
ensemble pq — e^^^^'^^Z, our proposal is to define an 
equivalent problem, such that the Hamiltonian H{q) is 
replaced by our functional -F(q) i. e. 

exp(-/3F(q)) 

= ^ (7) 

Here, as well as in standard statistical mechanics Z repre- 
sents the partition function ( this time Z = J e'^-'^^'^^dq) 
and P can be considered a parameter equivalent to the 
temperature. Automatically, this probability density 
maximizes a new "entropy" , and lead us to pose an equiv- 
alent state equation: 

A = F(q) - ST (8) 

Where A can be thought as the corresponding Helmholtz 
potential of the equivalent problem. In case of solving 
the eigenvalue problem in vibration theory, for example, 
F{q) turns to be Rayleigh's quotient i?(q) ^°i^^ . A briefly 
mathematical support is given in appendix |A] which as- 
sures the mathematical possibility of the probability den- 
sity proposed. 

Summarizing, since our problem of solving the ODE 
has been transformed into an equivalent thermal system, 
all the mathematical framework of Statistical Mechanics 
could be applied. For the sake of brevity, here we only 
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outline how several typical thermodynamical equilibrium 
variables such as entropy and specific heat, could be cal- 
culated using previous relations. For example, for the 
specific heat C(T) we have, by definition 



•^i^i- — — 



dT 



where () means equilibrium or ensemble average and 
could be calculated for each T by Monte-Carlo sampling 
(see section irvj) . 

Also the entropy could be obtained in the same form 
applying the well known relation 



dS_ _ C{T) 

It ~ T 

which integrated, using equation ([9]), give us 

T 



(10) 



S{T)^S{T^)- j 



i(£M!_(£M!)l,r (n) 



different configuration will differ from each other in an 
infinitesimal quantity which means, from the notion of 
a variation of a functional, that we reach a minimum 
q* = u* (solution of the ODE). By this way, the algo- 
rithm is supposed to be capable of avoiding relative min- 
ima, providing the absolute minimum of the functional. 
For the eigenvalue problem, this allows only to obtain 
the fundamental mode and lowest eigenvalue. However, 
higher modes which correspond to relative minima, could 
be obtained in principle if we impose restrictions to the 
functional, such as symmetry or number of nodes; these 
topics will be explored in future works. 



A. Solution generation 

Here, we address the problem of generating candidate 
solutions q' from an existing one q (solution generation 
scheme). To attempt this end, we will first make a brief 
survey of the calculus of variations^ . The variation of a 

functional F{u) — I f{u,u''-^\...,u''^^)dx which depends 



on the fmiction u{x) and its derivatives of order v, u^"^"^ 
is 



5F{u, 



,(1) ^ 



(13) 



IV. PROPOSED ALGORITHM 



where by definition 



Having established this new equivalent problem, we ap- 
ply the whole framework of equilibrium statistical analy- 
sis to obtain the extreme value of F{q) . Our main tool for 
solving problems through this paper is based on Metropo- 
lis Monte-Carlo algorithmic. Briefly, this algorithm en- 
sures a way to sample canonical ensemble proposing a 
transition probability Pq\q which satisfies detailed bal- 
ance. For our case, it takes the form 



Pq 



= cxp 



AF 



q',qj 



T 



(12) 



where T = 1//3, h = 1. 

With this transition probability, we construct an al- 
gorithm that can sample our n— dimensional configura- 
tional phase space Tin, whose limit as n — > cxd tends to 
Hilbert space Ti.A- Each equilibrium state, which can be 
reached by waiting sufficient time (Monte-Carlo steps), 
is a function of parameter T in which all possible con- 
figurations differ from each other in ±T. Obviously, this 
also means that for a given T > there exist many equi- 
librium configurations and only for T = it reduces to 
one. 

The proposed algorithm takes the previous idea and 
uses it in the same way as in standard simulated 
annealing^. We've just exposed that equilibrium states 
are governed by parameter T. So, if one diminishes its 
value in a way so as to go through successive equilib- 
rium states, then, the most representative points of the 
ensemble will be sampled. In the limit as T ^ 0, each 



e(u)5{u) 



dx 



5udx 



and Su is the variation of the function u{x) and is equal 
to eri{x), e — !■ 0; also ri(x) is an arbitrary function satis- 
fying sufficiently smooth conditions^. 

Our solution generation scheme starts generating a 
Markov chain in which the next iterate is built from the 
previous one making the change q' ^ q + Sq. The pro- 
posal consists in taking Jq in agreement with the master 
equation^. We consider R steps (fixed number of steps) 
for a given T. For the transition probability to have a 
value of 1/2 we set Ai^(q', q) ~ SF; then, 



^q',q = ^ = cxp ( - 



A^qVq) 

T 



exp 



e(q)'^q 

T 



(14) 

where equation (|13p has been used. Now, it is possible 
to calculate Sq = {dui{xi), ...,Sui{xi), ...Sun{xn)} con- 
sidering the integral in (jl3p as a sum over the n sites 
of the lattice. According to this, and considering every 
5ui{xi) identical and independent of cc;, we can calculate 
an estimation of the jump Su that satisfies ([T4|) 



5ui{xi) = 5u 



ln(l/2)T 



(15) 



Then, 5q = ...,5u]. This form of calculating 5q is 
one possible way of setting an estimation to the jump to 
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build the next iterate. Of course, other possibiHties are 
feasible to be implemented. However, it has the advan- 
tage of efficiently exploring configurational phase space in 
order to reach equilibrium states avoiding local minima 
(if there exists one). 



B. Algorithm implementation 

It is a relevant problem to determine the initial tem- 
perature To to start the algorithm. The one proposed 
by Kirkpatrick^ for initial temperature guesses in stan- 
dard simulated annealing technique (SA), proves to be 
simple and very effective. Obviously, other means are 
possible but we chose it for simplicity (seei^ and refer- 
ences therein). It consists of conducting a pilot survey of 
the solution space in which all increases in the objective 
function (the functional F(u) ) are accepted. The initial 
temperature is then calculated from once the initial 
transition probability -Pq',qo has been given ( 0.9 ~ 0.8 
commonly selected), according to 



To 



A/ 



(16) 



where A/ is an average increase in the objective func- 
tion, F{u). According to this formula it is then possible 
to initiate the algorithm. To trigger it for the first time, 
a high T is chosen. Then, compute A/ for the first 
hundred of attempts and apply (fT6|) . 

In the implementation of the algorithm, we decided 
to lower the temperature T half an order of magnitude 
between different equilibrium states. Then, we perform 
10^ runs to compute different equilibrium averages for a 
given temperature to ensure we reach equilibrium. 

Since it is impossible to reach absolute zero as a limit- 
ing value of the temperature T, a stopping criteria must 
be used as T ^ 0. This criteria could be condensed in the 
following way: the process of cooling of the system ends 
when, making a choice of a significative digit of the value 
of the functional being minimized, this digit averages to 
zero during the cooling procedure, while the other, more 
significative digits don't change. 



V. THE LANGEVIN APPROACH 

An alternative treatment of the problem is to present it 
within a stochastic framework. To this end, we propose 
to model the system with the Langevin equation. The 
use of Langevin equation to stochastically sample an ar- 
bitrary field is not a new ideai^. However, the authors 
ignore that the same approach have been attempted to 
find solutions of ODEs. Briefly, Langevin equation is a 
stochastic differential equation whose coefficients are ran- 
dom with given stochastic properties'^. It defines u{x,t) 



as a stochastic process provided that an initial condition 
is added, u{x,0). The developed method to solve ODEs 
consists of proposing a Langevin equation for the scalar 
field u{x) in the simulation time t, this is 



du{x, t) 
dt 



= -1 



1 6F[u] 
i{T) Su 



+ ri{x,t) 



(17) 



where ri{x, t) is a Gaussian thermal noise with zero mean 
and that satisfies 

< ri{x, t)r]{x', t') >= 2D{T)S{x ~ x')d{t ~ t') (18) 

There is a certain arbitrariness on the selection of D{T) 
and £,(T) but the product i{T)D{T) must satisfy Ein- 
stein relation, i.e. £,{T)D{T) — k^T to give the correct 
thermodynamical averages. It can be proved that the 
probability distribution P(u, t) approaches the equilib- 



rium distribution P, 



eq 



To obtain 



the minimum of F{u) and consequently the solution of 
the ODE, we must solve (fT7|) starting at a high tempera- 
ture To and then slowing it down gradually, which means 
passing through successive equilibria, until we reach a 
temperature near zero. Following the same reasoning as 
before we will obtain, by this way, the absolute minimum 
of the functional. Of course, the implementation of this 
cooling schedule is analogous to the previously developed 
procedure, so there is no need to go into detail again. 

In the end, to make possible a numerical solution of the 
field equation p7p . we work with the same discretization 
scheme as above. This time, Langevin equation for the 
scalar lattice field ui{xi) looks like 



dui{xi,t) 
di 



= -1 



1 SF[ui] 
aT) Suiixi) 



+ ri{xi,t) l^l..n (19) 



VI. NUMERICAL SIMULATIONS 

In this section we present several ODEs which were 
solved applying the previously developed methods. The 
cases under study are three 

• Case 1: Helmholtz equation. 

• Case 2: Fourth order beam equation. 

• Case 3: Fourth order beam deflection under static 
load. 

The proposed methods provide, for the eigenvalue 
problem (cases 1 and 2), the absolute minimum of the 
functional (lowest eigenvalue ) and the corresponding 
eigenfunction (fundamental mode). For case 3, this sit- 
uation no longer exists since there is only one minimum 
of the functional and this represents the unique solution 
of the ODE. 
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Case 1 



For the case of Helmholtz equation, u" (x) + \u{x) — 0, 
the corresponding functional F{u) reduces to 



of the considered configurations are presented in table (|T| 
where the number of algorithm iterations is also shown 
for time's estimation calculations. 



F{u) 



u{x)'^dx 



(20) 



This result is obtained from equation ^ after an easy 
manipulation. Fixed boundary conditions will be consid- 
ered in this case, this means u{0) — u{L) = 0. The func- 
tional is named Rayleigh quotient -R(q), as said before. 
Physically, it can be derived from the principle of conser- 
vation of energy since, if the total energy of the system 
remains constant, one can assure that Vmax = Tmax- 

For a vibrating system, let Tmax = w^T^^j^.; that is, 
^maa; the maximum kinetic energy of the system dur- 
ing a cycle of motion, with the square of the natural fre- 
quency, w^, factored out. Thus, from the total constant 
energy requirement, one can write 



(21) 



which is exactly (l20l) with oj^ — A. 

For the algorithm implementation, we start with a ran- 
dom function qo and then applies the solution generation 
routine and the transition probability Pq\q of (jl2p to- 
gether with the cooling schedule till the final temperature 
criterion is reached. For this case we have discretized the 
1 dimensional lattice in 16 sites, i. e. n = 16. However, to 
compute the functional values F{u), we have employed a 
cubic spline interpolation scheme using Matlab^ routines 
between u/'s values which smooths the representation of 
u{x) to calculate F{u). 

For the algorithm to perform the iterations over differ- 
ent configurations, new configurations are settled from 
old ones with a jump estimation Sq— {Su, Su} that is 
calculated from equation 03 



du 



ln(l/2)r||u|p 
2G{u) 



(22) 



where G(u) - 
function and 



J{u" + \u)dx represents Galerkin's error 
= ^u'dx. Here e(u) ^ 2G(m)/||w||2 
as can be seen by comparing equations and 
Obviously, u"{x) -f \u{x) = is only satisfied by the 
solution of the ODE u* . 

The initial random function qo can be observed in 
figure ll]) as initial config. In the same figure, inter- 
mediate candidate solutions which correspond to differ- 
ent equilibrium situations for different temperatures are 
shown. The final configuration was chosen following the 
final temperature stopping criterium. Numerical values 



1.5 



0.5 



-0.5 



final config 
intial config 



N part=16 



=0.43971 1 
interm config 




10 



15 



FIG. 1: Initial, intermediate and final configurations for the 
Helmholtz equation with fixed ends, showing the convergence 
of the algorithm to the minimum of the functional, Npart = 
16. 



config 


R(u) 




T 




N iter 


initial 


0.396959 


1 


10" 


-4 





interm 1 


0.083343 


1 


10~ 


-4 


1 


lo-^* 


interm 2 


0.048746 


1 


10" 


-4 


2 


10^ 


interm 3 


0.045338 


1 


10" 


-4 


3 


10^ 


final 16 


0.043871 


1 


10" 


-6 


1 


lO'' 



TABLE I: Rayleigh's quotient for the configurations of figure 
[T] with Npart = 16. Temperature and number of algorithm 
iterations are also shown. 



The selection of n = 16 or hereafter Npart = 16 was 
not arbitrary. A few number of particles at the beginning 
of the run makes the algorithm to run faster at large val- 
ues of the prescribed F{u) (we can call it the "energy" of 
the system). As this "energy" goes down, the algorithm 
starts to slow down its convergence rate due to the fi- 
nite number of particles. At this stage, we increment the 
number of particles as many times as it is necessary to 
get further significant reductions. This process has been 
named "stable resizing" due to the "resize" of the num- 
ber of considered particles. Figure ([2]) shows the results 
for Npart = 64. The exact solution (fundamental mode) 
appear in full line and the function that results of the 
minimization process is shown in cross line. The numeri- 
cal values of Rayleigh quotient are presented in table (|TT| . 
It can be noted that the difference between the exact and 
obtained R{u) is below 0.01% which demonstrates the 
convergence of the method to the exact solution. 
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5 10 15 



X 

FIG. 2: Final configuration x (cross) and exact solution — 
(full line) for the Helmholtz equation with Npart = 64. 



config 


R(u) 


T 




TV iter 


final 64 


0.0438650 1 


10" 


-9 


5- 10* 


exact 


0.0438649 









TABLE II: Rayleigh's quotient for the final configuration 
of figure [2] with Npart = 64. Here, the difi'erence between 
the exact solution and the solution provided by the algorithm 
(final configuration) differs only in 0.01%. 

B. Langevin approach to the Helmholtz equation 

For the case of the Helmholtz equation u" {x)^\u{x) — 
0, we have already shown that the corresponding func- 
tional is F{u) = (u' (x))"^ dx / J^^ u{x)'^dx which, in its 
discretized version is 

n n 

F{ui)=Y.{\/uifSxi/Y,uf6xi (23) 

1=1 1=1 

where Vu; = — ui)/5xi. With this form of F{ui) 

and applying an Euler updating scheme^, it can be sub- 
stituted in to give 



f(x) 




5 10 15 



FIG. 3: Initial, intermediate and final configurations for the 
Helmholtz equation with fixed ends for the simulations with 
the Langevin approach, showing the convergence of the algo- 
rithm to the minimum of the functional, Npart — 16. 



config 


R(u) 




T 




A'' iter 


initial 


0.273213 











interm 1 


0.115614 


1 


10- 


'2 


1 


10* 


interm 2 


0.053995 


1 


10" 


3 


1 


10* 


interm 3 


0.044929 


1 


10" 


4 


1 


10* 


interm 4 


0.044005 


1 


10" 


5 


1 


10* 


interm 5 


0.043875 


1 


10" 


6 


1 


10* 


final 16 


0.043866 


1 


10" 


7 


1 


10* 



TABLE III: Rayleigh's quotient for the configurations of fig- 
ure [3] with Npart = 16 for simulations with Langevin ap- 
proach. Temperature T and number of algorithm iterations 
are also shown. 



A stable resizing routine could also be applied in this 
case. The results for Npart = 64 are shown in table pvp . 



config 


R{u) 


T 


TV iter 


final 64 


0.0438661 1 


10-** 


1 ■ 10* 


exact 


0.0438649 







TABLE IV: Rayleigh's quotient for the final configuration of 
the simulations with Langevin approach with Npart = 64; the 
exact value of Rayleigh's quotient is shown for comparison. 



ui{t + St) = ui{t) - St 



where 



(ui+i - 2ui + ui-x)/Sxi + XuiSxi 



C. Case 2 

Tty/2D{T)7Ji{t) 

For the fourth order beam equation, u^^^^x) — Xu{x) 
0, the corresponding functional is 



^ufSxi 



1=1 



Fixed boundary conditions are also considered in this 
case. Numerical values of the obtained configurations are 
presented in table (jllip where the number of algorithm 



F(u) 



{u"ix))^dx 
u{x)'^dx 



(24) 



Clamped boundary conditions are considered for this 



iterations is also shown for time's estimation calculations. case, i. 



u(0) 



u{L) 



du{x)/dx\x=o 
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du{x)/dx\x=L = 0. As it was done before, the jump 
estimation can be calculated from equation (|15p result- 
ing in the same expression as (|22|) . The only change is in 
e{u) = 2G(u)/||u||2. This time, G{u) is 



G{u) ^ / - Xu)dx 



(25) 



We start the algorithm with a sinusoidal function, qo , 
which is not the solution of the differential equation in 
this case (see Figure S]). This was made to illustrate 
that the algorithm converges to the solution of the ODE 
independently of the initial condition that has been em- 
ployed. In figure ^ we also show the convergence of 
the algorithm to the solution of the beam equation as 
temperature T goes down. The number of particles is 
Npart = 16 for the first step of the process as was ex- 
plained before. Numerical values are shown in table (jVTl . 



config 


R(u) 




T 


N iter 


initial 


0.0372191 


5- 


10"'= 





interm 1 


0.0168120 


1 ■ 


10-^ 


2 • 10* 


interm 2 


0.0144376 


7.5 


•10-^ 


4-10* 


interm 3 


0.0123598 


2.5 


• 10"" 


6-10* 


final 16 


0.0098999 


1 • 


10-y 


2 • lO'' 



TABLE V: Rayleigh quotient for the configurations shown in 
figure |4] (fourth order beam equation); Npart — 16. Temper- 
ature T and the number of iterations are also shown. 



of =0.009887688 Npart =128 
exac 



A -0.5 




15 



FIG. 5: Final configuration x (cross) and exact solution 
(full line) of a beam with clamped ends; Npart = 128. 



D. Case 3 

For the fourth order beam deflection under a static 
load the equation is u^-^^^{x) — g{x). We have selected 
a concentrated load of the form g{x) ^ d{x ^ Xa) being 
S(x), the Dirac delta function. This time, the functional 
is 



03^=0.00989996 


Npart=16 


03^^^^=0.00988768 




interm config 




^/^^^v^ final config 




initial config ^^^Ny"^^ 









FIG. 4; Initial, final and intermediate configurations as the 
algorithm runs towards the solution of the beam equation 
with clamped ends. Npart = 16. 

To improve the precision of the algorithm and accel- 
erate its convergence, the number of particles were in- 
creased five times {Npart ~ 128) with respect to its ini- 
tial value. The final results can be observed in figure ([5]) 
where the Rayleigh quotient R is also shown. It is no- 
ticeable the better accuracy of the attained solution due 
to the "stable resizing" . 



F{u) = / {u"{x)fdx -2 u{x)g{x)dx (26) 



Simply supported boundary conditions are selected 
this time, u(0) = u{L) — d'^u{x)/dx^\x=o = 
d'^u{x)/dx'^\x=L = 0. In this case, the functional simply 
represents twice the potential energy in a beam subjected 
to the transverse force g{x). Applying equation (jl4p one 
is able to obtain the estimated jump which results 



5u 



\ii{l/2)T + 2du{xa) 
2G{u) 



(27) 



where G(u) = / u^^^^dx. This equation states that the 
estimated jump could be obtained once 5u{xa) is pro- 
vided. Since we have no knowledge of this quantity, it 
is reasonable to propose to be of order of Su. With this 
assumption, equation (|?7|) results 



5u 



ln(l/2)T 

'2{G{u) 



1 



(28) 



The algorithm was initialized with an initial configu- 
ration of the form of a sine function with random noise. 
This configuration, together with intermediates and final 
configurations are shown in figure Npart = 16 was 
chosen again since it provides good initial results. Nu- 
merical values of the different configurations are shown 
by table (jVip for Xa = 0.5L, where we rename I{u) = 
Fiu) . 



8 



config 


I{u) 




T 




Niter 


initial 


1.7231 


1 


10" 







interm 1 


0.1172 


1 


10" 


-4 


1 • 10^ 




u.uuiyyo 


1 
i 


1 0" 
lU 


-4 


• iU 


interm 3 


-0.0038 


1 


10" 


4 


1 ■ 10" 


interm 4 


-0.00808 


1 


10" 


-6 


2 • 10" 


interm 5 


-0.01103 


1 


10" 


-7 


1 • lO'^ 


final 


-0.011334 


1 


10" 


-8 


2 • 10"= 



TABLE VI: Functional value, /(it), for the case of the elastic 
beam equation under a static force of delta type with simply 
supported boundary conditions. 
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FIG. 6: Initial, final and intermediate configurations in the 
convergence to the solution of a beam under static load; 
Npart — 16. 



If a better precision wants to be reached, more particles 
must be included. Nevertheless with Npart = 16 the 
difference between the exact and final solution is below 
3%. 



VII. CONCLUSIONS 

A statistical approximation to the solution of ODEs 
were presented, in particular we solved some problems 
related with vibration theory. The developed method es- 
tablishes a way to attain the solution of an ODE proposed 
to solve physical problems, from statistical mechanics. 
The developed idea consists of thinking the domain of 
any functional (representing the ODE) as microstates of 
a thermal system (particle system) in contact with a heat 
bath at temperature T, that interacts trough a thermo- 
dynamic potential represented by the functional. Then, 
obtaining the minimum of the thermodynamic potential 
of this equivalent problem when T ^ is the same as to 
get the solution of the ODE. If only one solution exists, 
our method provides the unique solution of the ODE. 
In case we treat an eigenvalue equation, where infinite 



solutions exist, we obtain the absolute minimum of the 
corresponding functional or fundamental mode, as well 
as the fundamental eigenvalue. 

Three examples have been provided to show the con- 
vergence of method to the solution of the ODEs; the 
first two of them, represents a typical eigenvalue prob- 
lem, as it is the case for the Helmholtz equation and the 
fourth order beam equation, with infinite possible solu- 
tions. In these cases the method have shown to converge 
to the fundamental mode. For the third case, the beam 
equation under a static load, the method converge to 
its unique solution. Additionally, a resizing of the corre- 
sponding one dimensional lattice called "stable rezising" , 
which increases the discretization, was used to provide 
further reductions of the functional accelerating its con- 
vergence and improving its precision. 

In the same spirit, a general approach to the problem 
within an stochastic framework was presented with sim- 
ilar success. We proposed an equivalent Langevin equa- 
tion for the scalar field u in its way to find the minimum 
of the functional, solution of the ODE. Numerical exper- 
iments have been run that showed that the method also 
converges to the exact solution for the Helmholtz equa- 
tion. Of course, the other two cases, analyzed with the 
first developed method, can be treated within the same 
scheme but, for the sake of brevity, we postpone them for 
future works. 

Summarizing, the present article does not pretend to 
claim a place as an alternative method to other numeri- 
cal, probably faster and more precise, methods. Its virtue 
consists in that it is possible to establish a link between 
the problem of finding the solution of ODEs and an inter- 
acting particle system which can be treated by statistical 
mechanics; in this sense, we showed that some typical 
thermodynamical equilibrium variables, such as entropy 
and specific heat, could be in principle calculated. This 
provides a different perspective and new insights in solv- 
ing ODEs. 
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APPENDIX A 

When modelling a thermal system in equilibrium, the 
supportive mathematical structure is simple and not ex- 
clusive of physical systems. In this mathematical struc- 
ture, the entropy S is defined over a convex subset 
of R™"*"^, where a point of R™"''^ is denoted by 
{Xq, Xi....Xm) so it can be defined a function 
S": E^R 
such that^: 



• S' is concave 

• II > 0, being E = Xa 

• 5 is positively homogeneous of degree 1 

For a probabilistic model there is also a measure needed, 
then, the requirement is fulfilled if: 

• there is a class A defined by ^ = {p : ^ (0, oo)} 
( in physical systems, a point of U, is called a mi- 
crostate) 

• p is TT — measurable^ (tt is the reference measure) 

• p rfTT = 1 , ( p is the density of the microstate 
measure pdn) 



E{X, p) =< X >= /j^ Xp dn, {X is called an ob- 
servable) . 



There is a well known theorem.^ which ensures that 
exp • X) 



Z 



(Al) 



maximizes S{p) for all possible p £ A. Here, /3 G R™+^ 
and 



Z = / e-^*-^ dn 
In 



(A2) 
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